band <- individual_level %>% filter(abs(forcing) < 2500)
model_df1_1 <- lm(outcome1 ~ opium_legal + forcing + opium_legal*forcing + rugged + soil_quality + prec + temp, data=band)
model_df1_2 <- lm(outcome2 ~ opium_legal + forcing + opium_legal*forcing + rugged + soil_quality + prec + temp, data=band)
model_df1_3 <- lm(outcome3 ~ opium_legal + forcing + opium_legal*forcing + rugged + soil_quality + prec + temp, data=band)
model_df1_4 <- lm(outcome4 ~ opium_legal + forcing + opium_legal*forcing + rugged + soil_quality + prec + temp, data=band)

model_df1_1 <- coeftest(model_df1_1, vcov=cluster.vcov(model_df1_1, cluster = band$mkid14))
model_df1_2 <- coeftest(model_df1_2, vcov=cluster.vcov(model_df1_2, cluster = band$mkid14))
model_df1_3 <- coeftest(model_df1_3, vcov=cluster.vcov(model_df1_3, cluster = band$mkid14))
model_df1_4 <- coeftest(model_df1_4, vcov=cluster.vcov(model_df1_4, cluster = band$mkid14))

model_df1_1 <- tidy(model_df1_1) %>% mutate(outcome = "outcome1", level = "2.5km")
model_df1_2 <- tidy(model_df1_2) %>% mutate(outcome = "outcome2", level = "2.5km")
model_df1_3 <- tidy(model_df1_3) %>% mutate(outcome = "outcome3", level = "2.5km")
model_df1_4 <- tidy(model_df1_4) %>% mutate(outcome = "outcome4", level = "2.5km")

band <- individual_level %>% filter(abs(forcing) < 5000)
model_df2_1 <- lm(outcome1 ~ opium_legal + forcing + opium_legal*forcing + rugged + soil_quality + prec + temp, data=band)
model_df2_2 <- lm(outcome2 ~ opium_legal + forcing + opium_legal*forcing + rugged + soil_quality + prec + temp, data=band)
model_df2_3 <- lm(outcome3 ~ opium_legal + forcing + opium_legal*forcing + rugged + soil_quality + prec + temp, data=band)
model_df2_4 <- lm(outcome4 ~ opium_legal + forcing + opium_legal*forcing + rugged + soil_quality + prec + temp, data=band)

model_df2_1 <- coeftest(model_df2_1, vcov=cluster.vcov(model_df2_1, cluster = band$mkid14))
model_df2_2 <- coeftest(model_df2_2, vcov=cluster.vcov(model_df2_2, cluster = band$mkid14))
model_df2_3 <- coeftest(model_df2_3, vcov=cluster.vcov(model_df2_3, cluster = band$mkid14))
model_df2_4 <- coeftest(model_df2_4, vcov=cluster.vcov(model_df2_4, cluster = band$mkid14))

model_df2_1 <- tidy(model_df2_1) %>% mutate(outcome = "outcome1", level = "5km")
model_df2_2 <- tidy(model_df2_2) %>% mutate(outcome = "outcome2", level = "5km")
model_df2_3 <- tidy(model_df2_3) %>% mutate(outcome = "outcome3", level = "5km")
model_df2_4 <- tidy(model_df2_4) %>% mutate(outcome = "outcome4", level = "5km")

band <- individual_level %>% filter(abs(forcing) < 7500)
model_df3_1 <- lm(outcome1 ~ opium_legal + forcing + opium_legal*forcing + rugged + soil_quality + prec + temp, data=band)
model_df3_2 <- lm(outcome2 ~ opium_legal + forcing + opium_legal*forcing + rugged + soil_quality + prec + temp, data=band)
model_df3_3 <- lm(outcome3 ~ opium_legal + forcing + opium_legal*forcing + rugged + soil_quality + prec + temp, data=band)
model_df3_4 <- lm(outcome4 ~ opium_legal + forcing + opium_legal*forcing + rugged + soil_quality + prec + temp, data=band)

model_df3_1 <- coeftest(model_df3_1, vcov=cluster.vcov(model_df3_1, cluster = band$mkid14))
model_df3_2 <- coeftest(model_df3_2, vcov=cluster.vcov(model_df3_2, cluster = band$mkid14))
model_df3_3 <- coeftest(model_df3_3, vcov=cluster.vcov(model_df3_3, cluster = band$mkid14))
model_df3_4 <- coeftest(model_df3_4, vcov=cluster.vcov(model_df3_4, cluster = band$mkid14))

model_df3_1 <- tidy(model_df3_1) %>% mutate(outcome = "outcome1", level = "7.5km")
model_df3_2 <- tidy(model_df3_2) %>% mutate(outcome = "outcome2", level = "7.5km")
model_df3_3 <- tidy(model_df3_3) %>% mutate(outcome = "outcome3", level = "7.5km")
model_df3_4 <- tidy(model_df3_4) %>% mutate(outcome = "outcome4", level = "7.5km")


band <- individual_level %>% filter(abs(forcing) < 10000)
model_df4_1 <- lm(outcome1 ~ opium_legal + forcing + opium_legal*forcing + rugged + soil_quality + prec + temp, data=band)
model_df4_2 <- lm(outcome2 ~ opium_legal + forcing + opium_legal*forcing + rugged + soil_quality + prec + temp, data=band)
model_df4_3 <- lm(outcome3 ~ opium_legal + forcing + opium_legal*forcing + rugged + soil_quality + prec + temp, data=band)
model_df4_4 <- lm(outcome4 ~ opium_legal + forcing + opium_legal*forcing + rugged + soil_quality + prec + temp, data=band)

model_df4_1 <- coeftest(model_df4_1, vcov=cluster.vcov(model_df4_1, cluster = band$mkid14))
model_df4_2 <- coeftest(model_df4_2, vcov=cluster.vcov(model_df4_2, cluster = band$mkid14))
model_df4_3 <- coeftest(model_df4_3, vcov=cluster.vcov(model_df4_3, cluster = band$mkid14))
model_df4_4 <- coeftest(model_df4_4, vcov=cluster.vcov(model_df4_4, cluster = band$mkid14))

model_df4_1 <- tidy(model_df4_1) %>% mutate(outcome = "outcome1", level = "10km")
model_df4_2 <- tidy(model_df4_2) %>% mutate(outcome = "outcome2", level = "10km")
model_df4_3 <- tidy(model_df4_3) %>% mutate(outcome = "outcome3", level = "10km")
model_df4_4 <- tidy(model_df4_4) %>% mutate(outcome = "outcome4", level = "10km")

plot <-
  bind_rows(model_df1_1, model_df1_2, model_df1_3, model_df1_4,
            model_df2_1, model_df2_2, model_df2_3, model_df2_4,
            model_df3_1, model_df3_2, model_df3_3, model_df3_4,
            model_df4_1, model_df4_2, model_df4_3, model_df4_4) %>%
  mutate(outcome_label = case_when(outcome == "outcome1" ~ "Place of worship",
                                   outcome == "outcome2" ~ "Live in village",
                                   outcome == "outcome3" ~ "Live in neighborhood",
                                   outcome == "outcome4" ~ "Rent room",
                                   TRUE ~ NA_character_)) %>%
  mutate(level = factor(level, levels = c("10km", "7.5km", "5km", "2.5km"))) %>%
  filter(term == "opium_legal") %>%
  #mutate(Inequality = level) %>%
  ggplot(aes(x=estimate, y = outcome_label, group = level, shape = level, color = level)) +
  geom_point(position = position_dodge(width = 0.4)) +
  geom_errorbarh(aes(xmin=estimate-1.67*std.error, xmax = estimate+1.67*std.error), position = position_dodge(width = 0.4), height = 0) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  xlab("Effect of Opium Concession System") +
  theme_bw() +
  scale_color_grey() +
  theme(axis.ticks.y.left = element_blank(),
        panel.grid.minor = element_blank(), 
        axis.title.y = element_blank(),
        panel.grid.major.x = element_blank(),
        legend.title = element_blank(),
        axis.line.y.left = element_blank(),
        axis.line = element_line(colour = "black"),
        panel.border = element_blank(),
        panel.background = element_rect(fill = "transparent",colour = NA),
        plot.background = element_rect(fill = "transparent",colour = NA))


ggsave("./_4_outputs/figures/appendix_5_land_controls_plot.png", plot = plot,  width = 7, height = 4)
